Turbine Collision Risk: Overview Plots Distributions

Cédric Scherer https://cedricscherer.com (Leibniz Institute for Zoo and Wildlife Research)https://www.izw-berlin.de/en/home.html , Christian Voigt (Leibniz Institute for Zoo and Wildlife Research)https://www.izw-berlin.de/en/home.html
2022-02-03

Preparation

Show code
library(tidyverse)
library(here)
library(ggtext)
library(spatstat.core)
library(colorspace)
library(here)
library(patchwork)

theme_set(theme_void(base_size = 18, base_family = "Open Sans"))
theme_update(
  plot.title = element_markdown(hjust = .5, face = "bold"),
  plot.title.position = "plot",
  plot.margin = margin(rep(10, 4))
)

## source funcrions to generate distributions of bat passes
source(here("R", "src", "generate-distributions.R"))

Research Question

Is the predictive power of the current acoustic monitoring sufficient enough to estimate the true number of killed bats?

Predictions

The predictive power of the acoustic monitoring for the extrapolation of the expected number of stroke victims is decreasing

Open Questions

Inputs

We have the following parameters we can vary:

Show code
plot_dist <- function(data, title = NULL, color = NULL, shape = NULL, size = 1) {

  df_rotor <- as.data.frame(
    disc(radius = 1, centre = c(0, 0),  mask = FALSE, npoly = 5000)
  )
  
  area_monitored <- pi * unique(data$prop_monitored) * 2   ## pi *1^2 -> area 
  ## times 2 because we only use the lower half of the circle = half of the area later
  radius_monitored <- sqrt(area_monitored / pi)
  
  df_monitored <- 
    as.data.frame(
      disc(radius = radius_monitored, centre = c(0, 0), mask = FALSE,  npoly = 5000)
    ) %>% 
    filter(y <= 0)
  
  base <- 
    ggplot(df_rotor, aes(x, y)) + 
    geom_path(color = "black", alpha = .33, size = 2.5) + 
    coord_fixed() + 
    scale_x_continuous(limits = c(-1, 1)) +
    scale_y_continuous(limits = c(-1, 1))
  
  if (!is.null(color) & !is.null(shape)) {
    
     base + 
       geom_path(data = df_monitored, color = desaturate(lighten(color, .25), .25), size = 1.5) +
       geom_point(data = data, aes(color = monitored, shape = fatality, size = fatality, alpha = fatality, stroke = 1.2)) +
       scale_color_manual(values = c("grey20", color), guide = "none") + 
       scale_shape_manual(values = c(1, shape), guide = "none") +  
       scale_size_manual(values = c(2, 10), guide = "none") +  
       scale_alpha_manual(values = c(.3, 1), guide = "none") + 
       ggtitle(title)
    
  } else if (!is.null(color) & is.null(shape)){ 
    
    base + 
       geom_path(data = df_monitored, color = desaturate(lighten(color, .25), .25), size = 1.5) +
       geom_point(data = data, aes(color = monitored), alpha = .7, size = size) +
       scale_color_manual(values = c("grey20", color), guide = "none") + 
       ggtitle(title)
    
  } else {
    
    base + 
      geom_point(data = data, alpha = .7, size = size) + 
      ggtitle(title)
    
  }
}

Comparison Distribution Scenarios

Show code
m <- .2    ## proportion monitored
n <- 500L  ## number of basses
s <- 1L    ## seed

## uniform
df_uni  <- generate_distribution(distribution = "uniform", prop_monitored = m, n = n, seed = s, var = FALSE)

## random
df_ran  <- generate_distribution(distribution = "random",  prop_monitored = m, n = n, seed = s, var = FALSE)

## inner > outer (basic)
df_inn1  <- generate_distribution(distribution = "inner",   prop_monitored = m, n = n, seed = s, var = FALSE, skewness = 1)

## inner > outer (basic)
df_inn3  <- generate_distribution(distribution = "inner",   prop_monitored = m, n = n, seed = s, var = FALSE, skewness = 3)

## outer > inner (basic)
df_out1  <- generate_distribution(distribution = "outer",   prop_monitored = m, n = n, seed = s, var = FALSE, skewness = 1)

## outer > inner (basic)
df_out3  <- generate_distribution(distribution = "outer",   prop_monitored = m, n = n, seed = s, var = FALSE, skewness = 3)

## bottom > top (basic)
df_bot1 <- generate_distribution(distribution = "bottom",  prop_monitored = m, n = n, seed = s, var = FALSE, skewness = 1)

## bottom > top (extreme)
df_bot3 <- generate_distribution(distribution = "bottom",  prop_monitored = m, n = n, seed = s, var = FALSE, skewness = 3)

## top > bottom (basic)
df_top1 <- generate_distribution(distribution = "top",     prop_monitored = m, n = n, seed = s, var = FALSE, skewness = 1)

## top > bottom (extreme)
df_top3 <- generate_distribution(distribution = "top",     prop_monitored = m, n = n, seed = s, var = FALSE, skewness = 3)
Show code
gg_uni_r  <- plot_dist(df_uni,  title = "Uniform")
gg_ran_r  <- plot_dist(df_ran,  title = "Random")
gg_inn1_r <- plot_dist(df_inn1, title = "Inner > Outer")
gg_inn3_r <- plot_dist(df_inn3, title = "Inner > Outer (skewed)")
gg_out1_r <- plot_dist(df_out1, title = "Outer > Inner")
gg_out3_r <- plot_dist(df_out3, title = "Outer > Inner (skewed)")
gg_bot1_r <- plot_dist(df_bot1, title = "Bottom > Top")
gg_bot3_r <- plot_dist(df_bot3, title = "Bottom > Top (skewed)")
gg_top1_r <- plot_dist(df_top1, title = "Top > Bottom")
gg_top3_r <- plot_dist(df_top3, title = "Top > Bottom (skewed)")

## panel
p_r <- (gg_uni_r | gg_inn1_r | gg_out1_r | gg_bot1_r | gg_top1_r) / (gg_ran_r | gg_inn3_r | gg_out3_r | gg_bot3_r | gg_top3_r) +
  plot_annotation(tag_levels = "a", tag_suffix = ")")

p_r
Show code
ggsave(here("plots", "distributions", "distributions_all_rotor.pdf"), width = 24, height = 11, device = cairo_pdf)

## panel
p_r_l <- 
  wrap_plots(gg_uni_r, gg_ran_r, gg_inn1_r, gg_inn3_r, gg_out1_r, 
             gg_out3_r, gg_bot1_r, gg_bot3_r, gg_top1_r, gg_top3_r, ncol = 2) * 
  theme(plot.margin = margin(5, 40, 5, 40)) + 
  plot_annotation(tag_levels = "a", tag_suffix = ")")

ggsave(here("plots", "distributions", "distributions_all_rotor_long.pdf"), width = 12, height = 24, device = cairo_pdf)
Show code
gg_uni_m  <- plot_dist(df_uni,  color = "orange2", title = "Uniform")
gg_ran_m  <- plot_dist(df_ran,  color = "orange2", title = "Random")
gg_inn1_m <- plot_dist(df_inn1, color = "orange2", title = "Inner > Outer")
gg_inn3_m <- plot_dist(df_inn3, color = "orange2", title = "Inner > Outer (skewed)")
gg_out1_m <- plot_dist(df_out1, color = "orange2", title = "Outer > Inner")
gg_out3_m <- plot_dist(df_out3, color = "orange2", title = "Outer > Inne (skewed)r")
gg_bot1_m <- plot_dist(df_bot1, color = "orange2", title = "Bottom > Top")
gg_bot3_m <- plot_dist(df_bot3, color = "orange2", title = "Bottom > Top (skewed)")
gg_top1_m <- plot_dist(df_top1, color = "orange2", title = "Top > Bottom")
gg_top3_m <- plot_dist(df_top3, color = "orange2", title = "Top > Bottom (skewed)")

## panel
p_m <- (gg_uni_m | gg_inn1_m | gg_out1_m | gg_bot1_m | gg_top1_m) / (gg_ran_m | gg_inn3_m | gg_out3_m | gg_bot3_m | gg_top3_m) +
  plot_annotation(tag_levels = "a", tag_suffix = ")")

p_m
Show code
ggsave(here("plots", "distributions", "distributions_all_monitored.pdf"), width = 24, height = 11, device = cairo_pdf)
Show code
gg_uni_f  <- plot_dist(df_uni,  color = "orange2", shape = 4, title = "Uniform")
gg_ran_f  <- plot_dist(df_ran,  color = "orange2", shape = 4, title = "Random")
gg_inn1_f <- plot_dist(df_inn1, color = "orange2", shape = 4, title = "Inner > Outer")
gg_inn3_f <- plot_dist(df_inn3, color = "orange2", shape = 4, title = "Inner > Outer (skewed)")
gg_out1_f <- plot_dist(df_out1, color = "orange2", shape = 4, title = "Outer > Inner")
gg_out3_f <- plot_dist(df_out3, color = "orange2", shape = 4, title = "Outer > Inner (skewed)")
gg_bot1_f <- plot_dist(df_bot1, color = "orange2", shape = 4, title = "Bottom > Top")
gg_bot3_f <- plot_dist(df_bot3, color = "orange2", shape = 4, title = "Bottom > Top (skewed)")
gg_top1_f <- plot_dist(df_top1, color = "orange2", shape = 4, title = "Top > Bottom")
gg_top3_f <- plot_dist(df_top3, color = "orange2", shape = 4, title = "Top > Bottom (skewed)")

## panel
p_f <- (gg_uni_f | gg_inn1_f | gg_bot1_f | gg_top1_f) / (gg_ran_f | gg_inn3_f | gg_out3_f | gg_bot3_f | gg_top3_f) +
  plot_annotation(tag_levels = "a", tag_suffix = ")")

p_f
Show code
ggsave(here("plots", "distributions", "distributions_all_fatalities.pdf"), width = 24, height = 11, device = cairo_pdf)

Plots Example Distributions

Show code
seed <- 123L
distribution <- "inner"
color <- "orange2"

Small radius (30m) and low-frequency bats (20 kHz)

Show code
## small radius (30m) and low-frequency bats (20 kHz) 
##   -> small n             -> larger prop_monitored
small <- 100L
low_30 <- .5

df_30_low <- generate_distribution(n = small, prop_monitored = low_30, distribution = distribution, seed = seed, var = FALSE, skewness = 2)

(gg_d_30_low  <- plot_dist(df_30_low, size = 2, title = "30m radius + low-frequency"))
Show code
ggsave(here("plots", "distributions", paste0("distribution_example_", distribution, "_30m_low_freq_rotor.pdf")), width = 5, height = 5.2, device = cairo_pdf)

(gg_m_30_low <- plot_dist(df_30_low, size = 2, color = color, title = "30m radius + low-frequency"))
Show code
ggsave(here("plots", "distributions", paste0("distribution_example_", distribution, "_30m_low_freq_monitored.pdf")), width = 5, height = 5.2, device = cairo_pdf)

(gg_f_30_low <- plot_dist(df_30_low, color = color, shape = 4, title = "30m radius + low-frequency"))
Show code
ggsave(here("plots", "distributions", paste0("distribution_example_", distribution, "_30m_low_freq_fatalities.pdf")), width = 5, height = 5.2, device = cairo_pdf)

Small radius (30m) and high-frequency bats (40 kHz)

Show code
## small radius (30m) and high-frequency bats (40 kHz) 
##   -> small n             -> smaller prop_monitored
high_30 <- .17

df_30_high <- generate_distribution(n = small, prop_monitored = high_30, distribution = distribution, seed = seed, var = FALSE, skewness = 2)

(gg_d_30_high <- plot_dist(df_30_high, size = 2, title = "30m radius + high-frequency"))
Show code
ggsave(here("plots", "distributions", paste0("distribution_example_", distribution, "_30m_high_freq_rotor.pdf")), width = 5, height = 5.2, device = cairo_pdf)

(gg_m_30_high <- plot_dist(df_30_high, size = 2, color = color, title = "30m radius + high-frequency"))
Show code
ggsave(here("plots", "distributions", paste0("distribution_example_", distribution, "_30m_high_freq_monitored.pdf")), width = 5, height = 5.2, device = cairo_pdf)

(gg_f_30_high <- plot_dist(df_30_high, color = color, shape = 4, title = "30m radius + high-frequency"))
Show code
ggsave(here("plots", "distributions", paste0("distribution_example_", distribution, "_30m_high_freq_fatalities.pdf")), width = 5, height = 5.2, device = cairo_pdf)

Large radius (60m) and low-frequency bats (20 kHz)

Show code
## large radius (60m) and low-frequency bats (20 kHz)  
##   -> large n             -> larger prop_monitored
large <- 400L
low_60 <- .23

df_60_low <- generate_distribution(n = large, prop_monitored = low_60, distribution = distribution, seed = seed, var = FALSE, skewness = 2)

(gg_d_60_low <- plot_dist(df_60_low, size = 2, title = "60m radius + low-frequency"))
Show code
ggsave(here("plots", "distributions", paste0("distribution_example_", distribution, "_60m_low_freq_rotor.pdf")), width = 5, height = 5.2, device = cairo_pdf)

(gg_m_60_low <- plot_dist(df_60_low, size = 2, color = color, title = "60m radius + low-frequency"))
Show code
ggsave(here("plots", "distributions", paste0("distribution_example_", distribution, "_60m_low_freq_monitored.pdf")), width = 5, height = 5.2, device = cairo_pdf)

(gg_f_60_low <- plot_dist(df_60_low, color = color, shape = 4, title = "60m radius + low-frequency"))
Show code
ggsave(here("plots", "distributions", paste0("distribution_example_", distribution, "_60m_low_freq_fatalities.pdf")), width = 5, height = 5.2, device = cairo_pdf)

Large radius (60m) and high-frequency bats (40 kHz)

Show code
## large radius (60m) and high-frequency bats (40 kHz)  
##   -> large n             -> smaller prop_monitored
high_60 <- .04

df_60_high <- generate_distribution(n = large, prop_monitored = high_60, distribution = distribution, seed = seed, var = FALSE, skewness = 2)

(gg_d_60_high <- plot_dist(df_60_high, size = 2, title = "60m radius + high-frequency"))
Show code
ggsave(here("plots", "distributions", paste0("distribution_example_", distribution, "_60m_high_freq_rotor.pdf")), width = 5, height = 5.2, device = cairo_pdf)

(gg_m_60_high <- plot_dist(df_60_high, size = 2, color = color, title = "60m radius + high-frequency"))
Show code
ggsave(here("plots", "distributions", paste0("distribution_example_", distribution, "_60m_high_freq_monitored.pdf")), width = 5, height = 5.2, device = cairo_pdf)

(gg_f_60_high <- plot_dist(df_60_high, color = color, shape = 4, title = "60m radius + high-frequency"))
Show code
ggsave(here("plots", "distributions", paste0("distribution_example_", distribution, "_60m_high_freq_fatalities.pdf")), width = 5, height = 5.2, device = cairo_pdf)

Comparison Skewness

Show code
plot_skewed <- function(dist, skewness) {
  g <- plot_distribution(
    generate_distribution(distribution = dist, skewness = skewness, 
                          n = 300L, prop_monitored = .01, seed = 123L),
    title = paste0("Skewness = ", skewness), color = "black", print = FALSE
  )
  
  g <- g + theme(plot.title = element_markdown(size = 15, face = "plain"))
  
  return(g)
}

skewness_levels <- c(1:5, 10, 15, 20)
n <- length(skewness_levels)

## inner > outer 
plots_inner <- purrr::map2(rep("inner",  n), skewness_levels, ~plot_skewed(dist = .x, skewness = .y))  ## c(1, 3, 5, 7, 9, 15, 30, 60)
wrap_plots(plots_inner, ncol = 4) + plot_annotation(title = "Inner > Outer")
Show code
ggsave(here::here("plots", "distributions", "distributions_gradient_inner.pdf"), width = 12, height = 7, device = cairo_pdf)

## outer > inner
plots_outer <- purrr::map2(rep("outer",  n), skewness_levels, ~plot_skewed(dist = .x, skewness = .y))
wrap_plots(plots_outer, ncol = 4) + plot_annotation(title = "Outer > Inner")
Show code
ggsave(here::here("plots", "distributions", "distributions_gradient_outer.pdf"), width = 12, height = 7, device = cairo_pdf)

## top > bottom
plots_top <- purrr::map2(rep("top",  n), skewness_levels, ~plot_skewed(dist = .x, skewness = .y))
wrap_plots(plots_top, ncol = 4) + plot_annotation(title = "Top > Bottom")
Show code
ggsave(here::here("plots", "distributions", "distributions_gradient_top.pdf"), width = 12, height = 7, device = cairo_pdf)

## bottom > top
plots_bottom <- purrr::map2(rep("bottom", n), skewness_levels, ~plot_skewed(dist = .x, skewness = .y))
wrap_plots(plots_bottom, ncol = 4) + plot_annotation(title = "Bottom > Top")
Show code
ggsave(here::here("plots", "distributions", "distributions_gradient_bottom.pdf"), width = 12, height = 7, device = cairo_pdf)

Session Info
Show code
[1] "2022-02-03 11:53:59 CET"
Show code
R version 4.1.2 (2021-11-01)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19043)

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252   
[3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
[5] LC_TIME=German_Germany.1252    
system code page: 65001

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods  
[7] base     

other attached packages:
 [1] patchwork_1.1.1     colorspace_2.0-2    spatstat.core_2.3-2
 [4] rpart_4.1-15        nlme_3.1-153        spatstat.geom_2.3-1
 [7] spatstat.data_2.1-2 ggtext_0.1.1        here_1.0.1         
[10] forcats_0.5.1       stringr_1.4.0       dplyr_1.0.7        
[13] purrr_0.3.4         readr_2.0.2         tidyr_1.1.4        
[16] tibble_3.1.6        ggplot2_3.3.5       tidyverse_1.3.1    

loaded via a namespace (and not attached):
 [1] fs_1.5.0              spatstat.sparse_2.1-0 lubridate_1.8.0      
 [4] httr_1.4.2            rprojroot_2.0.2       tools_4.1.2          
 [7] backports_1.2.1       bslib_0.3.1           utf8_1.2.2           
[10] R6_2.5.1              mgcv_1.8-38           DBI_1.1.2            
[13] withr_2.4.3           tidyselect_1.1.1      downlit_0.4.0        
[16] compiler_4.1.2        textshaping_0.3.6     cli_3.1.0            
[19] rvest_1.0.2           xml2_1.3.2            labeling_0.4.2       
[22] sass_0.4.0            scales_1.1.1          askpass_1.1          
[25] goftest_1.2-3         systemfonts_1.0.3     digest_0.6.29        
[28] spatstat.utils_2.3-0  rmarkdown_2.11        pkgconfig_2.0.3      
[31] htmltools_0.5.2       highr_0.9             dbplyr_2.1.1         
[34] fastmap_1.1.0         rlang_0.4.12          readxl_1.3.1         
[37] rstudioapi_0.13       farver_2.1.0          jquerylib_0.1.4      
[40] generics_0.1.1        jsonlite_1.7.2        distill_1.3          
[43] magrittr_2.0.1        Matrix_1.3-4          Rcpp_1.0.7           
[46] munsell_0.5.0         fansi_0.5.0           abind_1.4-5          
[49] lifecycle_1.0.1       stringi_1.7.5         yaml_2.2.1           
[52] grid_4.1.2            crayon_1.4.2          deldir_1.0-6         
[55] lattice_0.20-45       splines_4.1.2         haven_2.4.3          
[58] tensor_1.5            gridtext_0.1.4        hms_1.1.1            
[61] knitr_1.36            pillar_1.6.4          markdown_1.1         
[64] reprex_2.0.1          glue_1.4.2            evaluate_0.14        
[67] pdftools_3.0.1        qpdf_1.1              modelr_0.1.8         
[70] vctrs_0.3.8           tzdb_0.1.2            cellranger_1.1.0     
[73] gtable_0.3.0          polyclip_1.10-0       assertthat_0.2.1     
[76] cachem_1.0.6          xfun_0.27             broom_0.7.11         
[79] ragg_1.1.3            memoise_2.0.1         ellipsis_0.3.2